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ABSTRACT 

We demonstrate that the formation of collapsing cores in subcritical clouds is 
accelerated by nonlinear flows, by performing three-dimensional non-ideal MHD 
simulations. An initial random supersonic (and trans-Alfvenic) turbulent-like 
flow is input into a self-gravitating gas layer that is threaded by a uniform mag- 
netic field (perpendicular to the layer) such that the initial mass-to-flux ratio is 
subcritical. Magnetic ambipolar diffusion occurs very rapidly initially due to the 
sharp gradients introduced by the turbulent flow. It subsequently occurs more 
slowly in the traditional near-quasistatic manner, but in regions of greater mean 
density than present in the initial state. The overall timescale for runaway growth 
of the first core(s) is several x 10 6 yr, even though previous studies have found a 
timescale of several x 10 7 yr when starting with linear perturbations and similar 
physical parameters. Large-scale supersonic flows exist in the cloud and provide 
an observationally testable distinguishing characteristic from core formation due 
to linear initial perturbations. However, the nonlinear flows have decayed suf- 
ficiently that the relative infall motions onto the first core are subsonic, as in 
the case of starting from linear initial perturbations. The ion infall motions are 
very similar to those of neutrals; however, they lag the neutral infall in directions 
perpendicular to the mean magnetic field direction and lead the neutral infall in 
the direction parallel to the mean magnetic field. 

Subject headings: ISM: clouds - ISM: magnetic fields - MHD - diffusion - tur- 
bulence 
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Introduction 



In the standard model of low-mass star formation, molecular clouds have an initially 
subcritical mass-to-flux ratio, and spend a relatively long time (~ 10 7 years) undergoing 
quasi-static fragmentation by ambipolar diffusion until cores come to be supercritical. After 
a core becomes supercritical, it collapses on a dyna mical timescale ( ~ 10 6 years) to form 
either one or a small multiple system of stars (e.g. 



Shu et al. 


1987; 


Mouschovias 


1991) 



Models of nonaxis ymmetric fragmentation by this mechanism have b e en presented recently 



in several papers (llndebetouw fc Zweibell l2000l ; iBasu fc Cioleki 12004 ; ICiolek fc Basul 12006 



Kudoh et al.l 120071 . hereafter KBOY). This scenario is supported by the observed relatively 
low efficiency of star formation in molecular clouds, i.e. , only 1 — 5% for clouds as a whole but 



several times greater within cluster- forming cores (see lLada fc Ladal 120031 ) . An alternative 
model is that star formation begins in clo uds with supercritical mass-to -flux ratio that have 
additional support due to turbulence (e.g., Nakancfl998 ; Hartmann 2001 ). In this model, star 



formation occurs relatively rapidly, as turbulence dissipates over a dynamical timescale, and 
the large-scale magnetic field plays a minor role. The rapid star formation model is supported 
by the observational results that the age spreads of young stars in nearby molecular clouds 
are often a few x 10 6 years (lErmegreenll2000l ; lHartmannll200ll . 120031 ) . and the low fraction of 
clouds that are observed to be in a pre-star-formation state. 



Recently, iLi fc Nakamural (12004 ) and lNakamura fc Lil (120051 ) have modeled mildly sub- 
critical clouds and shown that the timescale of cloud fragmentation is reduced by supersonic 
turbulence. They showed this by performing two-dimensional simulations in the thin-disk 
approximation. Such a model can explain both relatively rapid star formation and the rel- 
atively low star formation efficiency in molecular clouds that is not well explained if star 
formation starts from a supercrit i cal cl oud. In this Letter, we study the three-dimensional 
extension of the ILi fc Nakamural (120041 ) model, by including a self-consistent calculation of 
the vertical structure and dynamics of the cloud. We confirm that supersonic flows can 
significantly shorten the timescale of core formation in three-dimensional subcritical clouds, 
and clarify how the scaling of various physical quantities allows this to occur. 



2. Numerical Model 



We solve the three-dimensional magnetohydrodynamic (MHD) equations including self- 
gravity and ambipolar diffusion, assuming that neutrals are much more numerous than ions. 
Instead of sol ving a detailed energy equa tion, we assume isothermality for each Lagrangian 
fluid particle kudoh fc Basull2003l I2OO6}) For the neutral-ion collision time and associated 
quantities, we follow IBasu fc Mouschoviasl (I1994I ). The basic equations are summarized in 
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KBOY. 

As an initial condition, we assume hydrostatic equilibrium of a self-gravitating one- 
dimensional cloud along the z-direction in a Cartesian coordinate system (x,y,z). Though 
nearly isothermal, a molecular cloud is usually surrounded by warm material, such as neutral 
hydrogen gas. Hence, we assume that the initial sound speed makes a transition from a low 
value Cgo to a high value c sc at a distance z = z c , with a transition length z& (see eq. 16 of 
KBOY). We take c 2 c = 10c^ , z c = 2H , and z& = 0.1H , where c s0 is the initial sound speed 
at z = 0, H = c s0 /y/2nGpo, and po is the initial density at z — 0. A numerical sol ution for 



the initial density distribution shows that it is almost the same as the lSpitzerl (119421 ) solution 
for an equilibrium isothermal layer in the region < z < z c . We also assume that the initial 
magnetic field is uniform along the z-direction. 

A set of fundamental units for this problem are c s o, H , and po- These yield a time 
unit to = Hq/c s o. The initial magnetic field strength (B ) introduces a dimensionless free 
parameter 

_ 8vrp _ 8np c 2 s0 _ c 2 s0 

the ratio of gas to magnetic pressure at z — 0. In the above relation, we have also used 
Vao = B / y/Anpo, the initial Alfven speed at z — 0. In the sheet-like equilibrium cloud with 
a vertical magnetic field, Po is related to the mass-to-flux ratio for Spitzer's self-gravitating 
layer. The mass-to-flux ratio normalized to the critical value is ps = 2%G 1 ^ 2 Tis/Bo, where 
£s = 2p H is the column density of the Spitzer cloud. Therefore, Po = p|. Although the 
initial cloud we used is not exactly the same as the Spitzer cloud, Po is a good indicator 
of whether or not the magnetic field can prevent gravitational instability. For the model 
presented in this Letter, we choose p = 0.25 so that the cloud is slightly subcritical. For 
this choice, Vao — 2.8 c s0 . Dimensional values of all quantities can be found through a choice 
of po and c s o- For example, for c s o = 0.2 km s -1 and n = po/rrin = 10 4 cm -3 , we get 
H = 0.05 pc, t = 2.5 x 10 5 yr, and B = 40 pG if p = 0.25. 

The level of magnetic coupling in the partially ionized gas is characterized by numerical 
values of the ion number density rij and neutral-ion collision timescale r n i. From eqs. (8) and 
(9) of KBOY, and using standard values of parameters in that paper as well as the values of 
units used above, we find an initial midplane ionization fraction x^o = h^o/^o = 9.5 x 10~ 8 
and a corresponding neutral-ion collision time r n i,o = 0.1 lto- The ionization fraction x\ and 
timescale r n j at oth er densities can b e found from the initial midplane values given that they 



both scale oc p 1 ^ 2 (lElmegreenl 1 1 9791 ) . 



In this equilibrium sheet-like gas layer, we input a nonlinear perturbation to v x and 
at each grid point. Independent realizations of v x and v y are generated in Fourier space 
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with amplitudes drawn from a Gaussian distribution and consistent with power spectrum 
v\ oc k~ A . Here, k = (k 2 + k 2 ) 1 ^ 2 is the total wavenumber. Our adopted power spectrum 
means that the large scale has mor e energy than the sma ll scale. This perturbation is 
consistent with that implemented in iLi fc Nakamural (12004] ) . The rms value of the initial 
velocity perturbation in physical space, v a , is about 3c so , so that v a — Vao as well. 



Kudoh et al. 


1999; 


Oeata et al. 


2004) 



Fl, \y 



< z < AHq, with a number of grid points for each direction (N x , N y , N z 



BOY (see also 
| < 8ttH and 
(64,64,40). 



3. Results 



The top panel of Figure 1 shows the time evolution of the maximum density p max at 
z = 0. The simulation is stopped when p max = 30p . The solid line shows the result 
when we input an initially nonlinear supersonic perturbation to the (3q = 0.25 model. The 
dashed line ((3q = 0.25) and the dash-dotted line (@q = 4) show the result in the case of the 
small (linear) initial perturbation that we have applied in KBOY. This figure shows that the 
timescale of collapsing core formation for the nonlinear perturbation case is much shorter 
than that for the linear perturbation case, when /3 is the same. Even when the initial 
cloud is subcritical (/% = 0.25), the core formation occurs on almost the same timescale as 
that of the supercritical (/3q = 4) linear perturbation case 0. The bottom panel of Figure 
1 shows the time evolution of the maximum value of density at z = (p m ax) and /3 at the 
location of maximum density (/3 pmax )- At first, /3 pmax increases rapidly up to ~ 0.9 due to 
rapid ambipolar diffusion in the highly compressed regions caused by the initial supersonic 
perturbation. However, there is enough stored magnetic energy in the compressed region 
that it rebounds and starts oscillations, with /3 pmax around 0.7 and increasing gradually. 
Eventually, /? pmax becomes > 1 and the dense region collapses to form a core. This figure 
implies that j3 is a good indicator to see whether a subregion of the cloud is supercritical 
or not. The evolution of p max confirms that there is an initial compression followed by a 
rebound to a lower density (still greater than the initial background value) and subsequent 
oscillations until a runaway collapse starts when continuing ambipolar diffusion has created 
a region with (3 > 1. 

Figure 2 shows an image of the logarithmic density at the last snapshot (t = 20.5to, 
when p max = 30po) for the nonlinear perturbation case of the subcritical cloud ({3q = 0.25). 



^^When the initial cloud is supercritical and the perturbation is supersonic, the collapsing core formation 
happens quickly, at t ~ to, from the initial flow. This may be too rapid to agree with observations. 
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The top panel shows the cross section at z — 0, and the bottom panel shows the cross 
section at y — —5.9H . The value of y for the bottom panel is chosen so that the vertical 
cut passes through the maximum density point. A collapsing core is located in the vicinity 
of x — —20H ,y = —5H . The size of the core is similar to that created by linear initial 
perturbations (see Fig. 2 in KBOY), although the shape is notably less circular. 

The top panel of Figure 3 shows the density and x-velocity of neutrals and ions along 
an x-axis cut at y — —5.9H , z — 0, taken from the snapshot illustrated in Figure 2. The 
x-velocities show infall motion toward the center of the core, although the core itself is 
moving with nonzero negative x-velocity. The relative infall speed to the core is subsonic 
and about 0.35c s o- It is comparable to the case of initial linear perturbation (see Fig. 9 in 
KBOY). However, there are systematic motions throughout the simulation region that are 
still supersonic at this time, even though the initial turbulent energy has decayed somewhat. 
This is qualitatively distinct from the corresponding case with linear perturbations, and 
should be observationally testable. For the single core that is formed in this simulation, the 
systematic core motion is about 0.5c s o in the x-direction. However, it can be even larger 
in other realizations. The ion velocity Vi can be related to the neutral velocity v and the 
magnetic acceleration a M in the assumed limit of low ion inertia by the relation 

Vi = v + a M r ni . (2) 

Figure 3 reveals that the ion infall motions are smaller in magnitude than the neutral motions 
due to the retarding magnetic acceleration; however the relative drift in the x-direction 
between ions and neutrals is typically within 0.05c s o- In order to see the force balance in the 
core, we plot the x-accelerations in the bottom panel of Figure 3 along the same x-axis as the 
top panel. These values are normalized by the computational units. The inward gravitational 
acceleration (a G ) dominates around the core. The acceleration from the magnetic force (om) 
and thermal pressure force (or) resist the contraction. The net acceleration (a nct ) is working 
toward the center of the core. 

The top panel of Figure 4 shows the density and z-velocities along a ,2-axis cut at 
x = —20H ,y = — 5.9H from the snapshot illustrated in Figure 2. The z-velocities also 
show infall motion toward the center of the core. The relative neutral infall speed to the 
core is also subsonic and about — 0.35c s o- The ion velocity in the vertical direction is almost 
identical to the neutral velocity but actually very slightly greater in magnitude at the peak. 
The bottom panel of Figure 4 shows the ^-accelerations along the same z-axis as in the top 
panel. In the z-direction, a G is nearly balanced by a T , and a M is almost negligible, but 
working in the same direction as gravity. The net acceleration (a ne t) is also working toward 
the center of the core. The inward pointing vertical magnetic acceleration associated with 
an hourglass-type magnetic field morphology explains why the ion velocity in the z-direction 
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is greater in magnitude than the neutral velocity. 



Discussion and Conclusions 



Our three-dimensional MHD simulations have shown that the supersonic nonlinear flows 
significantly reduce the timescale of collapsing core formation in subcritical clouds. It is of 
order several x 10 6 years for typ ical parameters , or_~ 10 times less than found in the lin- 



ear initial perturbation studies of iBasu fc Ciolekl (120041 ). ICiolek fc Basul (120061 ). and KBOY. 
Our supersonic perturbations a re also tran s- Alfvenic, in agreement with analysis of ob- 
se rved magnetic field st r ength s (jBasul 120001 ) . It is also in line with the theoretical result 
of iKudoh fc Basul (120031 . 120061 ) that turbulent clouds will settle into a trans- Alfvenic state 
when global motion/expansion is allowed, even if driven with turbulence that is initially 
super-Alfvenic. 

To see how accelerated ambipolar diffusion can occur, we note that the magnetic in- 
duction equation (see, e.g., eq. 3 of KBOY) and our assumptions of ionization balance can 
be used to estimate the diffusion time oc p 3 ^ 2 L 2 /B 2 , where L is the gradient length scale 
introduced by the initial turbulent compression and B is the magnetic field strength. Be- 
cause the compression by the nonlinear flow is nearly one-dimensional, the magnetic field 
scales roughly as B oc L^ 1 within the flux freezing approximation. If the compression is 
rapid enough that vertical hydrostatic equilibrium cannot be established (unlike in previous 
calculations using the thin-disk approximation), then p oc L~ l as well (i.e., one-dimensional 
contraction without vertical settling), and oc L 5 / 2 . This means that diffusion can occur 
quickly (and lead to a rapidly rising value of (3) if the turbulent compression creates small val- 
ues of L. If diffusion is so effective during the first turbulent compression that a dense region 
becomes magnetically supercritical, then it will evolve directly into collapse. Alternately, 
the stored magnetic energy of the compressed (and still subcritical) region may lead to a re- 
expansion of the dense region. The timescale for this, in the flux-freezing limit, is the Alfven 
time ta, oc Lp 1 ! 2 /B, which scales oc L 3 / 2 for the above conditions. Thus, decreases more 
rapidly than ta, and sufficiently small turbulent-generated values of L may lead to enough 



magne tic diffusion that collapse occurs before any reexpansion can occur. See lElmegreen 



(120071 ) for some similar discussion along these lines. Ultimately, whether or not reexpansion 
of the first compression can occur depends on the strength of the turbulent compression, 
mass-to-flux ratio of the initial cloud, and neutral-ion collision time. In a preliminary study 
of such effects, we found that keeping v a and (3q fixed but allowing significantly poorer or 
stronger neutral-ion coupling yielded differing results. For r ni)0 = 0.3to (poorer coupling), a 
collapsing core formed due to the first compression, at t ~ 1.4t , while for r nii0 = 0.05to it 



-7- 



happened after more oscillations than in our standard model, at t ~ 85to- A full parameter 
study will be required to elucidate. 

If reexpansion of the initial compression does occur, as in the standard model presented 
in this paper, then there is enough time for the vertical structure to settle back to near- 
hydrostatic equilibrium, in which case B oc p 1 / 2 . Since the compressed and reexpanded region 
executes oscillations about a new mean density, it is convenient to analyze the scalings in 
terms of the density p. The diffusion time now scales as oc p -1 / 2 . This yields a scaling of 
that is the traditionally used one (and is satisfied by design in the thin-disk approximation). 
However, the diffusion occurs more rapidly than it would in the initial state due to the 
elevated value of p in the compressed but oscillating region (see Figure 1). 

SB was supported by a grant from NSERC. Numerical computations were done mainly 
on the VPP5000 at the National Astronomical Observatory. 
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Fig. 1. — Top: the time evolution of maximum densities at z — 0. The solid line shows the 
evolution for an initially nonlinear supersonic perturbation and (3q = 0.25. The dashed line 
(/3q = 0.25) and the dash-dotted line (/3q = 4) show the evolution for models with a linear 
initial perturbation as calculated by KBOY. The dotted line shows the evolution for an 
initially nonlinear supersonic perturbation and [3o = 0.25, but without ambipolar diffusion. 
Bottom: evolution of the maximum density (solid line) at z = of the simulation box for the 
model with fi Q = 0.25 and nonlinear initial perturbation, and evolution of j3 at the location 
of maximum density (dashed line). 
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Fig. 2. — Logarithmic density image at t = 20.5io for the nonlinear perturbation case of 
the subcritical cloud (/3 = 0.25). The top panel shows the cross section at z — 0, and the 
bottom panel shows the x — z cross section at y = —5.9H . An mpeg animation of the 
evolution of the density image from t = to t — 20.5to is available online. 
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Fig. 3. — Top: the density (solid line), x- velocity of neutrals (dashed line), and a:- velocity 
of ions (dotted line) along an x-axis cut at y = —5.9H , z = 0, in the snapshot shown in Fig. 
2. The ^-positions are measured by offset from x c = —20Hq, which is the maximum density 
point for the core. Bottom: the x- accelerations along the same x-axis as the top panel. The 
solid line is the gravitational acceleration (ac), the dotted line is the magnetic acceleration 
(om), the dashed line is acceleration due to thermal pressure (ay), and the dash-dotted line 
is the net acceleration (a net ). All quantities are normalized by the computational units. 



- 12 - 



100 F 




-/o r , , , i 

-2-1 1 2 



z 

Fig. A.— Top: the density and z-velocity of neutrals and ions along a z-axis cut at x — 
—20H ,y = —5.9H , in the snapshot shown in Fig. 2. Bottom: the ^-accelerations along the 
same z-axis as the top panel. The line styles are the same as those in Fig. 3. 



